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Elastic netwok models coarse grain proteins into a network of residue beads connected by springs. We 
add dissipative dynamics to this mechanical system by applying overdamped Langevin equations of motion 
to normal-mode vibrations of the network. In addition, the network is made heterogeneous and softened 
at the protein surface by accounting for hydration of the ionized residues. Solvation changes the network 
Hessian in two ways. Diagonal solvation terms soften the spring constants and ofF-diagonal dipolc-dipole 
terms correlate displacements of the ionized residues. The model is used to formulate the response functions 
of the electrostatic potential and electric field appearing in theories of redox reactions and spectroscopy. We 
also formulate the dielectric response of the protein and find that solvation of the surface ionized residues 
leads to a slow relaxation peak in the dielectric loss spectrum, about two orders of magnitude slower than the 
main peak of protein relaxation. Finally, the solvated network is used to formulate the allosteric response of 
the protein to ion binding. The global thermodynamics of ion binding is not strongly affected by the network 
solvation, but it dramatically enhances conformational changes in response to placing a charge at the active 
site of the protein. 

Keywords: Protein solvation, elastic network model, dielectric spectroscopy, redox reactions, allostery, dissi- 
pative dynamics 



I. INTRODUCTION 

Folding a globular protein in water largely places po- 
lar/ionized residues to its surface, while moving the non- 
polar residues to its core. The resulting structure is 
not unique, and a number of conformations with close 
energy minima always exist. Conformational changes 
are required for function. They are achieved by either 
populating the existing (quasi)stable states (sampling 
of pre-existing equilibria^— ) or by shifting the existing 
minimum-energy conformation to a new configuration 
minimum upon perturbation, such as ligand binding (in- 
duced fit mechanism^). 

Conformational transitions involve several types of free 
energy penalty. The free energy of elastic deformation 
relative to the native structure, involving global shape 
alteration of the protein, is the most prominent penaltyi^ 
This is not the only long-ranged component of the overall 
protein's thermodynamics since electrostatic interactions 
are also involved in several ways. Changing the protein 
conformation alters the interactions between its atomic 
charges, but also, to a significant extent, the free energy 
of solvation of these charges by hydration water. Water 
clearly affects the flexibility of proteins As a fast, highly 
polar subsystem, it follows adiabatically the large-scale 
protein motions, continuously stretching and loosening 
the protein structure by strong protein-water solvation 
forces. It lowers the barriers of transitions between the lo- 
cal minima of the rugged landscape at the energy bottom 
of the native basin of attraction,— accelerating the rate 
of conformational changes.— When dried, proteins stiffen 
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and their relaxation time, as probed by dielectric spec- 
troscopy, increases by about six orders of magnitude.-^ 

The goal of this paper is to develop an efficient com- 
putational algorithm to include hydration in calculations 
of conformational flexibility of large protein complexes. 
Our starting point is to coarse-grain the protein into an 
elastic network of beads, a formalism known as elastic 
network model (ENM).— These types of models aim 
at calculating the elastic energy of deformation near the 
equilibrium structure and the directions of normal-mode 
displacements corresponding to the slowest normal-mode 
vibrations. The typical coarse-graining is achieved on the 
scale of a single residue by replacing it with a single rigid 
bead. The beads are then connected by elastic springs 
physically capturing the connectivity, shape, and packing 
of the residues in the folded protein structure. 

The ENM coarse-graining of the elastic energy has 
proven to be very successful.—""— Global elastic defor- 
mations of a protein are mostly affected by its shape 
and mass distributioni^^ri^i Electrostatics is another good 
candidate for coarse-graining. Coulomb forces are long- 
ranged and effectively average out the variations of the 
local structure. The flnal outcome for the free energy 
of electrostatic interactions is mostly determined by the 
overall density and distribution of the protein charge and 
the dipolar polarization of the hydration water. This 
physical reality is addressed by generalized Born solva- 
tion models designing fast computational algorithms to 
calculate the free energy of electrostatic solvation!^ 

The problem addressed here is two-fold. First, we want 
to re-normalize the elastic network by water's hydra- 
tion. Given the large free energy of hydration of the 
surface residues, the network of beads is expected to 
be softer at the interface. We achieve this goal here 
by integrating out the dipolar polarization of the hy- 
dration water using formalisms developed in the liquid- 
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FIG. 1. Cartoon displaying the electrostatic perturbation in- 
duced by a half redox reaction of transferring electron to the 
heme in the protein used as an example. The elastic deforma- 
tion of the protein shifts the ionized surface residues, shown 
by charges at the surface, but also results in adiabatic move- 
ments of water dipoles solvating them (shown by arrows). 



state theory of polar liquids . i The result is an analyt- 
ical model, a solvated dissipative electro-elastic network 
model (sDENM), which assigns lower force constants to 
springs attached to ionized interfacial residues. The sec- 
ond issue is the calculation of the response functions re- 
lated to problems affected by the protein electrostatics. 
Here, we consider three types of problems: (i) electro- 
static response to a probe charge or dipole placed inside 
the protein, (ii) response of the protein to a uniform ex- 
ternal field (dielectric spectroscopy), and (iii) elastic re- 
sponse at a given site of the protein to altering the charge 
state of a distant residue (allosteric action). 

The first type of problems appear in redox reactions in- 
volving proteins^ and in optical and IR spectroscopy^Si^ 
when the position of a spectral line is affected by the local 
electric field. In redox reactions, electron is transferred 
by tunneling from an electron donor to an active site 
(shown as protein heme in Fig. [T]) . The dynamics of this 
processes, and the activation barrier required to produce 
resonant conditions for electron tunneling, can be calcu- 
lated from the electrostatic response function X4>{^)- It 
arises from an elastic deformation, shifting the atomic 
charges of the protein, caused by transferring the elec- 
tron, but also from the change in the positions and ori- 
entations of water dipoles hydrating the surface residues 
(Fig. [1]). For spectroscopic applications, it is the dipole 
moment of the chromophorc that is altered by light ab- 
sorption. The corresponding response functions xe(j-^) is 
the one of the electric field acting on the chromophorc 
dipole and responsible for spectral solvatochromismj^ 

The second problem addressed here is the dynamic 
susceptibility of the protein to a uniform external elec- 
tric field produced in the dielectric spectroscopy exper- 
iment. Dielectric spectroscopy of partially wet protein 
powders has identified a number of generic relaxation 
peaks, the assignment of which has been problematic— 
The solvent-renormalized network model developed here 
results in three relaxation peaks of the protein assigned 
to fast backbone vibrations (fastest), the global shape- 
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FIG. 2. Cartoon showing propagation of piezoelectric per- 
turbation caused by binding an ion to residue i and produc- 
ing a displacement of residue j. The electric force exerted 
by the ion is propagated throughout the protein as an elas- 
tic deformation indicated by chains of arrows. In contrast, 
the displacements of ionized surface residues are propagated 
as water-mediated, dipole-correlated surface motions. Ion- 
ized surface residues combine into a global, correlated net for 
transmitting signals, which does not necessarily require a spe- 
cific binding site. 



altering movements (main peak), and the slow motions 
of highly solvated charged residues with significant extent 
of solvent exposure (slowest). 

Finally, the last type of problems considered here is 
the allosteric response J^i^ It specifies the alteration 
in the structure of the protein produced by a perturba- 
tion at a distant sitcj^ The perturbation can be achieved 
by localized ligand binding, often carrying a charge 
(Fig. [5]). Allosteric signaling usually involves oligomeric 
proteins, although single-domain proteins also display 
allostery.iii^i^ Given that two equilibrium conformations 
are involved, two equilibrium sets of atomic coordinates 
need to be considered for a full description of allosteric 
signaling. The barrier to the transition between the two 
equilibrium structures is the free energy of the protein 
elastic deformation, which can be approximated as cross- 
ing of two harmonic elastic wellsi^i^"— The elastic free 
energy is quadratic as a function of global normal-mode 
displacements, but can change its functional form to a 
linear function when localized (cracking) excitations, cor- 
responding to local unfolding events, are producedi^ 

The lowest elastic barrier is reached along the lowest 
curvature path on the free energy surface against defor- 
mation, i.e., the lowest frequency of the elastic vibra- 
tion. The allosteric pathways are therefore often asso- 
ciated with the lowest frequencies of harmonic motions 
near the two equilibrium structures.— The dissipative dy- 
namics of these motions can therefore be explored in the 
framework of response functions referring to a single equi- 
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librium conformationi^ The formulation of such response 
functions in the framework of a dissipative electro-clastic 
network is our purpose here. 

We calculate the dynamics of displacement of a distant 
residue in response to changing the charge at the binding 
site of an allosteric protein. The main question here is 
how water modifies the response. Wc find that solvation 
of ionized residues provides a potential mechanism alter- 
native to the typically anticipated clastic propagation of 
the perturbation. 

Despite differences in packing and connectivity of 
residues in different regions of a folded protein, elastic 
response tends to be non-specific, spreading out over the 
entire volume of the protein (Fig. [2]). In contrast, a net 
of ionized residues potentially provides an alternative, 
surface-bound propagation of the perturbation by water- 
mediated allostery. The water-mediated cross-coupling 
between the displacements of ionized residue scales as 
r~^. Therefore, an alteration of the charge at a bindig 
site can propagate large distances over the network of 
surface residues, instead of, or an addition to, the bulk 
elastic deformation. The allosteric response can then be 
channeled to a site where a conformational change is re- 
quired for function. 



II. MODEL 

The Hamiltonian of the protein hydrated by polar wa- 
ter can be generally written in the following form 

i/ = i?(R)-^E,,.m,. (1) 

Here, i?(R.) is the solvent-unperturbed Hamiltonian of 
the protein depending on the manifold of atomic co- 
ordinates R. Further, E^j is the electric field acting 
from residue i of the protein on dipole moment nij of 
water. Both the protein coordinates R and the water 
dipoles raj fluctuate with the instantaneous configura- 
tion of the protein-water system; summation over all 
residues i = 1 , . . . , and all waters j = 1 , . . . , iV^ is 
taken in Eq. 

Several approximations need to be made in the tran- 
sition from the general Hamiltonian in Eq. ([1]) to a sol- 
vatcd elastic network. The first approximation is the 
assumption that an equilibrium configuration of the pro- 
tein atomic coordinates is available from the structural 
data, and quadratic expansion in atomic displacements 
can be done around it. The model thus deals with one 
conformational state of the protein only and has nothing 
to say about transitions between distinct protein confor- 
mations. 

The first step in coarse-graining the model is to re- 
place the collection of protein atomic coordinates with 
a collection of beads. We will follow here the standard 
approach^r^i of representing each residue with a single 
bead, with its position given by the coordinates of the 



C" atom. The quadratic expansion of i?(R) in small dis- 
placements Srf = rf — of individual beads relative 
to equilibrium positions ro,i leads to the relation 

E^iC/2)J2Hr/6rT6r^ (2) 

in which Hff is a 3N x 3N Hessian matrix and C is 
the scaling force constant; a, (3 indicate the Cartesian 
projections, and summation over repeated Greek indices 
is assumed. 

The electrostatic component of the problem is repre- 
sented by the standard formulation of atomic force fields. 
This implies that each atom of the protein carries the 
charge qik, where i = 1, . . . , N numbers the residues and 
k represents an atom within residue i. The linear, in the 
displacements Svi, expansion of the protein- water inter- 
action term in Eq. ([1]) results in the following equation 

6E, • nij = ^ QikSvi ■ • mj, (3) 

k 

where Tij ~ — ViVjjroi — r^l^^ is the dipolar tensor 
connecting the C" of residue i with the dipole of water 
j- 

In the elastic network constructed here all atoms of a 
given residue experience a uniform displacement Svi from 
their equilibrium positions. The librations of the residues 
arc therefore neglected. This approximation leads to a 
significant simplification in Eq. ([3]) since only charged 
residues with '^^.Qik = li ^ contribute to the sum. 
Clearly, uniform displacements of only charged residues 
contribute to the creation of the dipole moment fluctua- 
tion 6fii ~ QiSri. We therefore get for the energy of the 
protein-water system 

H = iC/2) HrfSrfSr^ - qM^^nfm]. (4) 

We now proceed to calculating the free energy of the 
hydrated protein by tracing out the fluctuations of the 
dipole moments of water. Adiabatic approximation is 
assumed at this step, that is we assume that water is a 
fast subsystem, equilibrating to each instantaneous con- 
figuration of the network of beads. Since only global, 
relatively slow motions of the protein are modeled by the 
elastic network, this approximation is expected to be ac- 
curate. 

Averaging over the configurations of the dipole mo- 
ments of water produces partial free energy, i.e., free 
energy depending on the manifold of instantaneous dis- 
placements bvi. If the fluctuations of the dipolar polar- 
ization field of water are Gaussian, this free energy is 
given by the following equation 

F^{Cl2)Y^Ht^hr'i8rl (5) 

where the Hessian matrix, renormalized by solvation, be- 
comes 

mt = ii^'-C~'^mr (6) 
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In Eq. §1, the rank-2 tensor representing the dipo- 

lar response of the solvent to a dipole J/Xj created at the 
residue j. The dipolar polarization field created in the 
solvent in response to this perturbation then propagates 
to induce the dipole SfJbi at residue i. The correspond- 
ing free energy cost contributes to the renormalization of 
the Hookean force constants for the residues involved, as 
represented by the second term in Eq. ([S]). 

The physical meaning of the solvation terms in Eqs. ([5]) 
and © is quite clear. At i = j, the second term in Eq. 
(|6|) represents the solvation free energy of the fluctuation 
dipole 5fii. Correspondingly, i ^ j terms are the water- 
mediated couplings of the dipolar fluctuations through 
the solvent polarization. Note that direct Coulomb in- 
teractions between Sfii and Sfij are not included in Kij. 
These electrostatic terms propagate the elastic perturba- 
tion through the solvent, by hopping between the ionized 
residues, in addition to the direct propagation of elastic 
forces through elastic contacts of neighbors in the net- 
work (Fig. [5]). We note that the standard coarse-grained 
models of protein electrostatics, such as generalized Born 
modclS)^ do not include off-diagonal terms in their sol- 
vation free energy. These terms are however sufficiently 
long-ranged, scaling as with the distance between the 
residues, and they can potentially modify the response of 
the hydrated protein to either mechanical or electrostatic 
perturbation. 



A. Polar response of hydration water 

The dipolar susceptibility Kij in Eq. ([6]) generally re- 
quires either liquid-state models of solvation or electro- 
static continuum approaches for its calculation. Here, we 
start with the former to introduce a sequence of steps to 
reduce the full complexity of polar response to a clear 
physical picture and a computationally efficient algo- 
rithm. 

The susceptibility Kij is given by the convolution of 
the dipolar tensors Ti = T(roi — r), representing the 
electric field of the dipole Sfii at the point r in water, 
with the spacial correlation function of the dipolar fluc- 
tuations of water interfacing the protein. It can be conve- 
niently represented by the convolution of inverted space 
k-integrals2^i2iiii 



/^,, =f,(k)*x(k,k')*T,(k') 



(7) 



Here, the asterisks between tensors represent tensor con- 
traction over common indices and integration over com- 
mon k-variables. Further, the response function x(k, k') 
depends on two wave-vectors to reflect the inhomoge- 
neous nature of the problem caused by the presence of 
the protein in solution. Finally, Ti(k) is the Fourier 
transform of the dipolar tensor taken over the volume 
n occupied by water 



T,(k) = / T(r-roO0(r)e 
In 



ik-r 



dr. 



(8) 



As we have shown elsewhere^^ the nonlocal part of 
x(k, k') is mostly due to transverse polarization fluc- 
tuations, given by the component of the dipolar polar- 
ization perpendicular to the unit vector k = k/fc i^^i"^^ 
This transverse response is in fact fairly small for most 
solvation problems^ (e.g., the Born solvation energy is 
entirely longitudinal) and will be neglected here. This 
approximation eliminates the dependence on the second 
wave-vector with the result^ 



X(k,k') = kkxhfc)(27r)3<5(k-k') 



(9) 



Here, xf^(fc) is the longitudinal dipolar susceptibility of 
the homogeneous liquid depending on the scalar magni- 
tude k only. It is typically given as a product of the 
density of dipolcs in the liquid y and the longitudinal 
structure factor S^{kyM^ xH^) = (Sy /4:TT)S^{k). The 
dipolar density parameter y = (At: / 9) (3 pm? is defined 
by the liquid number density p and molecular dipole m; 
(3 = \/{k^T) is the inverse temperature. 

With the form of the response function given by Eq. 
([9]), the convolution in Eq. ([7]) is reduced to a 3D in- 
tegral. While this problem is numerically tractablcj^i^ 
the number of integrals to be evaluated is ^ /2, where 
Ni is the number of ionized residues. This is still a nu- 
merically intense computation, and simplifications are 
desired. 

We will further simplify the problem by modeling the 
calculation of the dipolar tensor of a given residue in 
Eq. ([8]). The full calculation of the Fourier transform re- 
quires numerical integration over the volume outside the 
typically complex shape of the protein^Sl To avoid this 
computationally extensive step, the concept of relative 
accessible surface area^ will be employed here. Specifi- 
cally, the volume integral in Eq. ([8]) will be replaced with 
the integral outside the sphere of radius s, representing 
the average distance of the closest approach of the water 
molecules to the residue, and scaled with the fraction of 
the surface area ai exposed to the solvent 



T,(k) = -47rDa,:^iMe''^-"-°-. 



(10) 



In this equation, D = 3kk — 1, jn{x) is the spherical 
Bessel function, and ai is the ratio of the solvent exposed 
area a,- to the overall surface area of the residue 



a,/(47rs2). 



(11) 



The reduction of Eq. (|T0|) yields an analytical solution 
for the solvent response function 



a{s, r) 



pOij 



(12) 



Here, indices are dropped for brevity in r = and 

is the direct-space dipolar tensor connecting beads i and 
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j. The first summand in the brackets in Eq. ([T2|) repre- 
sents the solvation free energy of the dipole at a charged 
bead and the second term represents the dipole-dipolc 
interaction between the two charged beads. 

The coefficients a{s,r) and b{s,r) in Eq. (|12|) arc ob- 
tained as one-dimensional integrals including the longi- 
tudinal structure factor of the liquid (k) to account for 
non-local correlations between the solvent dipoles. These 
integrls are listed and calculated in the Appendix. We 
show there that the dependence on s and r can be lifted 
in these functions and they in fact are well represented 
by constants, a(s, r) = S^{0)A, b{s,r) = S^{0). Here, 
5^(0) = {3y)~^{l — ej^) represents the longitudinal di- 
electric response of a homogeneous polar liquid with the 
dielectric constant e^. We finally get for the solvation 
tensor 



1 



1 



A 



SapSij - (1 - Sij)T"'^ 



(13) 

The constant A = 3.54 is calculated in the Appendix 
assuming s = 4.4 A for the closest-approach distance be- 
tween the center of a surface amino acid and the oxygen 
of water. The result is not strongly affected by the choice 
of s. Note, however, that A absorbs the thermodynamic 
state of the solvent into it and will change with its ther- 
modynamic state (temperature, pressure, etc) through 
the corresponding alterations of the polarization struc- 
ture factor. All calculations and MD simulations pre- 
sented here refer to the temperature of 300 K and ambi- 
ent pressure. 



B. Dissipative elastic network 

The free energy of the hydrated protein in Eq. ([5]) is 
a quadratic form in residues' displacements Jr^. It can 
be used to calculate the response to external perturba- 
tions or equilibrium variances once the network Hessian 
H"j^ has been specified. We will use here the Hookean 
springs Hamiltonian suggested by Tirion<^ This poten- 
tial, E(R) = {C/2)J2tj Dij{rij - ro,ij)^, describes the 
elongation Vij = I'^i — rj\ between nodes i and j in the 
network characterized by one universal force constant C 
(''o,ii = |ro.i ^ '"ojD- In this potential, Dij is the connec- 
tivity matrix. Its value is set to unity when is within 
the cutoff distance Tc and is set to zero otherwise. In 
addition, Dij ~ e > 1 for covalently bound neighbors. 
This scaling accounts for a stronger bonding of residues 
along the backbone, and is known to better model the 
vibrational density of states of the protein4i Finally, the 
rcnormalization of the network by solvation of and elec- 
trostatic interactions between ionized residues will follow 
Eqs. dSD and (O. 

The equations of motion for the network beads need 
to be specified in order to calculate the time-dependent 
response functions4S. The elastic network obviously lacks 
dissipative dynamics typical for soft condense phases. Al- 



ternatives to purely mechanical equations of motions can 
be sought in terms of Langevin dynamics of individual 
beads fi^i^ Introducing dissipation at the level of indi- 
vidual beads is not necessarily an obvious choice?^ and 
we have previously opted to introduce dissipation to nor- 
mal modes qm diagonalizing the network Hessiaui^ The 
equation of motion for such overdampcd dynamics is^ 



C{t - t')q^{t')dt' + A„q,„ = F{t) + R(i), (14) 



where C,{t — t') is a memory function, F(t) ~ Ft^e*"* 
is an external oscillating force, and R(t) is a randomly 
fluctuating force. The latter satisfies the generalized 
fiuctuation-dissipation relations^^i^ 



(R(t))=0, (R(t)-R(O)) = fcBTC(<). 



(15) 



The eigenvalues Am of normal modes qm in this equa- 
tion arc obtained by diagonalizing the Hessian in Eq. 
©. They are therefore affected by solvation softening 
the interface. Wc indeed observe a shift of the vibra- 
tional density of states to softer modes when the network 
is solvated. 

Applying Laplace-Fourier transformi^ to Eq. re- 
sults in the displacement response function for the col- 
lective mode q^. It is given as a scalar function con- 
necting the average displacement to the external field^^ 
(q™(a;)) = Xm(w)Fc^, where (qm(t)) = (qm(a;))e"^*. 
From Eq. (fT4|) . one gets 



Xm(w) = iujC,{Lo) -f A„ 



(16) 



where C{u)) is the Laplace-Fourier transforms of the fric- 
tion kernel C,{t). The entire set of 3A^ eigenvalues Am is 
produced by diagonalizing the Hessian with the unitary 
matrix U. The inclusion of all eigenvalues of the Hessian 
results in the response function of the bead displacements 



a/3 



7/3 



(17) 



The distance- and self-correlation functions of pro- 
tein residues typically show two characteristic relaxation 
times of overdampcd motion and, correspondingly, two 
Debye peaks in their loss spectra. Therefore, following 
the prescription of our previous work;^ we use a two- 
Debye form of Xm{^), which two characteristic friction 
coefficients, C/ and Qh 



Xm(^^) 



1 - a 



iujQh + Am i^Ci + ^r. 



(18) 



where the amplitude a specifies the relative weight of 
each relaxation component. 



C. Electrostatic response functions 

The network response function Xij(aj) in Eq. (|17p de- 
scribes the displacement of residue i induced by a weak 
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oscillating force applied to residue j. Since the residue 
displacement uniformly moves all of it atomic charges, 
this linear susceptibility can be used to build electro- 
static response functions of either electrostatic potential 
or electric field at a given location within the proteini^ 
Assume that an oscillatory probe charge qo{t) ~ Qt^e"^* 
is placed at some location Tq within the protein. This 
charge will act on residue i with the force — go(^)Eoi, 
where Eq^ is the electric field produced at rg by all 
charges of residue i at their equilibrium positions 



E, 



O'i 



k 



(19) 



Here, qik are the atomic charges of residue i with the 
equilibrium coordinates r^fc. 

The force acting on residue i will propagate through 
the elastic network to residue j according to the response 
function Xiji^)- The displacement of that residue will in 
turn produce an alteration of the electrostatic potential 
of the protein at Tq. After summing over all residues 
in the network, one arrives at the frequency-dependent 
susceptibility of the electrostatic potential 



(20) 
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FIG. 3. Free energy surfaces representing the free energy 
penalty (reversible work) of changing the charge at a bindig 
site of an allosteric protein. The calculations are done for 
attaching carbamoylphospate (go2 ~ —2) to the bacterial 
enhancer-binding protein NtrC (Fig. [S]). The curves refer 
to DENM in the unphosphorylated state {qoi ~ 0) and to 
DENM and sDENM in the phosphorylated state {qo2 = -2). 
The elastic network is defined with ksT/C = 0.75 and 
e = 125; the cutoff radius is 15 A. The free energy of binding 
AF is unknown and was set at —0.5 eV for the purpose of 
illustration. 



susceptibility, one arrives at the relation 

SM'^{u) = J2q.q,X:fHE^^. (22) 



This susceptibility determines the alteration of the elec- 
trostatic potential produced by the charges of the protein 
matrix (50o(w) at the position Tq of the probe charge qi^: 
5(j)oluj) = X0(w)'7w 

Similarly, one can define the electric field alteration 
(5Eo produced by the protein matrix at the point Tq in 
response to placing an oscillating probe dipole /x^ at that 
point. This frequency-dependent susceptibility is based 
on convoluting (w) with the dipolar tensors T^^ ^ 
T(ro — Vik) connecting the residue charge qik, located at 
r^fc, to the position of the probe dipole at rg. The result 
is 



E 

i,j.k.l 



Ta-y 7(5 



xlj{^)T-l'qjl. 



(21) 



Here, as above, summation runs over the repeated Greek 
indices denoting Cartesian projections of the correspond- 
ing tensors. The difference in signs in Eqs. ((20)) and ((2T|) 
comes from the fact that the free energy invested into 
the creation of the potential alteration is (l/2)gtj(S0o(a;), 
while for the dipole one has — (l/2)/.t„ ■ 5Eq[ijj). 



D. Dielectric response 

When a uniform oscillatory external field Eo(t) = 
E^e*"* is applied to a protein, it induces the dipole mo- 
ment (5M(aj) = 5iJ,j{uj). Since only movements of 
the charged residues produce non-zero dipoles, 5Ts/1{lj) = 
'^jqj5vj[u). Substituting the network displacements 



Assuming that the external field is along the z-axis of the 
laboratory frame, one gets for the dipolar susceptibility 

Xm(w) - ^g.g,x-/(^) = (l/3)^g.g,xr(^)- (23) 



This dipolar susceptibility refers to the dipole moment 
induced at a single protein molecule. It can be used to 
calculate the complex- valued dielectric constant ep(a;) of 
a protein sample (either powder or polycrystal) by apply- 
ing the standard derivation of the theory of dielectricsi^ 
The result is 



{ep{uj) - l)i2e(Lo) + 1) _ 47r 
9ej,(a;) 3 



PpXm(w), (24) 



where pp is the number density of the protein molecules 
in the material. 



E. Allosteric response 

As an example of the application of the formalism of 
response functions to the allosteric response of a protein, 
we will consider the displacement 5rf{uj) of residue i in 
response to binding a charge q{t) = Qi^e"^* at position 
ro. The response is therefore effectively of piezoelectric 
typC)^ creating deformation at a distant site in response 
to electric stimulus. 

Binding of an ion causes both global and local per- 
turbations of the protein. From the global perspec- 
tive, it changes the free energy of the entire protein by 
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the free energy of binding AF and, in addition, exerts 
Coulomb forces acting on all charges of the protein. The 
global perspective can be studied, as is typically done 
in electrostatics;^ by asking what is the free energy cost 
Fi (q) of transferring a small probe charge q to the binding 
site, where i = 1,2 labels the two stable conformations 
of the protein, ion-free and ion-bound. The electrostatic 
potential susceptibility [Eq. (|20|) ] addresses this question. 
The free energy cost is obviously 



1 



F^iq)^-xf{0){q-q^or + F,O, 



(25) 



where gio = 0, Fio = and 920 = Qi is the charge of the 
binding ion and F20 = AF is the binding free energy. We 
have also added the dependence of the response function 
on the protein state since susceptibility ^ (0) can be 
sensitive to structural changes of the protein. 

The amount of transferred charge q can be viewed as 
the reaction coordinate for the global free energy cost of 
binding the partial charge q. The crossing of the free 
energy surfaces, Fi{q'^) = ^2(9^), will define the transi- 
tion state and the corresponding free energy barrier. This 
picture allows for a "Marcus inverted region behavior" ^ 
i.e. there is an optimal binding free energy minimizing the 
free energy barrier. The activation barrier starts to grow 
when AF falls below the optimal value (inverted region). 
In addition, because the curvatures of the two parabo- 
las may differ, there is a scenario in which no crossing 
in the inverted region occurs, i.e. the activation barrier 
becomes infinite and the reaction is not allowed. Note, 
however, that the curvatures of two surfaces calculated 
for the NtrC protein studied below are nearly identical 
(Fig. [3]). There is also little sensitivity of the overall 
free energy functions to solvation of the surface residues 
(compare DENM and sDENM calculations). This lack of 
global sensitivity is in stark contrast with a strong effect 
of solvation on individual residue displacements, as we 
show below. 

Most of the interest in the field is driven not by the 
global thermodynamics of binding, but by the need to 
understand biological function caused by iti^i^ A typi- 
cal problem is to calculate the displacement of a distant 
residue in response to binding. We will approach this 
question, as above, by considering an oscillatory charge 
placed at rg. This charge will interact with each residue j 
by the electric field Eqj given by Eq. (fTQ]) . That interac- 
tion creates the force acting on each bead in the network, 
SF^"{io) = -Ef^q^. 

In the linear response approximation^ the average 
displacement of residue i is given by summing up the 
forces produced by the charge q^i at all residues of the 
network with their response function xt!^ (w) propagat- 
ing the force at j into a displacement at i 



(26) 




MD 

DENM 

sDENM 

charged residues 
from B-factors 



50 

residue No. 



100 



FIG. 4. Mean-square displacements of cytochrome B562 
(cytB) C"'s from MD simulations, DENM, and sDENM cal- 
culations. The network parameters in the DENM/sDENM 
calculations are: ksT/C — 0.75 A^, e = 125, cut-off radius 
is 15 A. The solvent accessible surface for the loop residues 
labeled in Fig. [5] is scaled down to 45 A^. 



placement 



1/2 



(27) 



where 5r'f{uj) is the real part of the complex- valued dis- 
placement and summation over repeated Greek indices is 
performed. 

The frequency uj of the oscillatory charge might effec- 
tively represent the time-scale of charge binding, such as 
the frequency of binding/unbinding events, typically oc- 
curring on the nanosecond time-scale for small electrolyte 
ions.— The limit a; = in this formalism will represent 
stationary, i.e., adiabatically slow binding. 



III. RESULTS 

Before presenting the results of specific calculations, 
we start with some crude estimates of the effect of solva- 
tion of charged residues on the properties of the clastic 
network. Our previous calculations of electrostatic prop- 
erties of redox proteins were done with the elastic spring 
constant of C = 0.6 kcal/ (mol A^), consistent with other 
estimates in the literaturej^ii^ Given this force constant, 
one can estimate the effect of solvation on the network 
Hessian. We consider the diagonal element in Eq. (|15p , 



(28) 



In the calculations below we will present the scalar dis- 



With q = e, s = 4A A, = 78, and A = 3.54, the 
second term becomes 30a|. This estimate suggests that 
any singly-charged residue exposed to water to more than 
a quatcr of its surface will have a negative elastic con- 
stant with its non-covalent neighbors and will be held in 
the equilibrium position only by covalent bonds within 
the network. Such ionized residue would lose mechani- 
cal stability and dissolve in water if not held in place by 
its covalent neighbors. It is clear that solvation makes a 
major effect on the elastic response of charged residues. 
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FIG. 5. Cartoon of cytochrome B562 (cytB) showing the 
positions of C" (spheres) colored by residue charge: charged 
(green) and uncharged (red). The charged residues marked 
green (also green points in Fig. 3)) are also required to have 
greater than 0.16, used as a threshold number. The remaining 
C" are marked red. The side chain atoms are colored by 
charge: negative (red), positive (blue), and neutral (white). 
The heme iron is colored brown while the remaining atoms 
of the heme are blue. Numbers label the unstable residues in 
the loop for which the water-exposed surface was scaled down 
to 45 in order to maintain the network stability. 
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FIG. 6. Loss spectra Xb('^)/x'b(0) and X0(^)/x;(O) for cytB. 
Compared are MD, DENM, and sDENM calculations. To 
show the sensitivity of sDENM calculations to solvation of 
the loop residues in cytB, the results of choosing the solvent- 
accesible area of ai — 45 and of ai — 50 A^ are shown. The 
two-Debye relaxation parameters are — 30 ns, ("h — 0.006(i, 
and a = 0.35 [Eq. (|18p ]. The elastic network is defined with 
ksT/C = 0.75 A^ e = 125, and the cutoff radius of 15 A. 



A. Residue displacements and electrostatic 
response 

The standard approach of experimental verification 
and parameterization of elastic protein networks is to 
compare the root-mean-squarc displacements (rmsd's) of 
residues with experiment or Molecular Dynamics (MD) 
simulations. Crystallographic B-factors arc often used^ 
but those are of limited valueJ^ It was noted that re- 
ported B-factors are dominated by rigid-body motions 
of the proteins in the crystal)^ In addition, there is a 
clear mismatch between the reported rmsd's of proteins 
in crystals and in their flexibility in solution, as is illus- 
trated in Fig. [4] comparing rmsd's from B-factor of C"'s 
with their rmsd's found from MD. The MD simulations 
were done for hydrated cytochrome B562 (cytB, PDB 
entry 256B, Fig. [5]) as described elsewhere.— 

The mismatch between both the B-factors and the 
standard ENM as compared to MD is particularly no- 
table for the flexible loop (residues 46 to 55) contain- 
ing several ionized residues (Figs. [4] and [5|). This region 
is clearly not restricted to a single configuration in so- 
lution and instead wanders through a number of semi- 
stable conformations. The network, required to reside in 
a single conformation, is expected to lose stability be- 
cause of this and similar segments. The standard ENM 



clearly avoids this instability by over-restricting the fiexi- 
ble residues. In contrast, when renormalization by solva- 
tion is introduced in sDENM, the network loses stability, 
as expected, due to solvation of the loop residues labeled 
in Fig. [5l Since calculations cannot be performed with 
an unstable network, we have artificially restricted the 
network by scaling down the solvent-accessible area of 
the residues labeled in Fig. [5] from the values calculated 
with VMDffi (in the range 120-150 A^) to 45 A^. This 
rescaling prevents mechanical instability of the network, 
but preserves the physical reality of a flexible loop, as is 
seen from the corresponding rmsd's in Fig. [D 

Figure [6] shows the results of the calculations (cytB) for 
the electrostatic potential susceptibility x<t>i^) ^^'^ ^^e 
tensor contraction X£;(<^) — X's'i^) ^'^^ ^^'^ electric field 
susceptibility. The effect of solvating surface residues is 
less pronounced for these susceptibilities, in particular 
for the more long-ranged electrostatic potential. A slow 
relaxation component, not resolved on the length of the 
MD trajectory, appears for the electric field susceptibil- 
ity. This slow component arises from much slower mo- 
tions of highly solvated residues in the network, also seen 
in the dielectric response of the protein. 
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FIG. 7. xX/('^)/Xm(0) for cytB. The parameters of the net- 
work are the same as in Fig. [6l 
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B. Dielectric susceptibility of the protein 

The calculated imaginary part of the dielectric suscep- 
tibility (loss function) x'iii'^) in shown in Fig. [71 Simi- 
larly to the case of x'^(ix)), it clearly shows the emergence 
of a slow peak, about two orders of magnitude slower 
than the main peak. The slow component comes from 
the hydrated residues with high exposure to water. This 
is clearly seen from the sensitivity of the slow peak to 
the solvent-accessible surface assigned to the residues of 
the loop. We need to note that the network completely 
neglects librations of polar residues, focusing only on the 
polarization fluctuations produced by translational mo- 
tions of the charged residues. An additional dielectric in- 
tensity might therefore come from the components miss- 
ing from the model. 

Experimentally, partially hydrated protein powders 
show three relaxation processes at low temperatures, 
which merge into two processes at ambient temperaturcj^ 
The fastest and the slowest processes disappear when the 
protein is dried. The relaxation time of the main peak 
from dielectric measurements matches well the relaxation 
time from neutron scattering experiments, in which the 
protein and hydration water signals can be separated 
by deuteration. The main peak is therefore assigned to 
global protein motionsi^ In this regard, the main peak in 
Fig. [7] can be tentatively aligned with the main peak of 
dielectric measurements. 

The slowest observed peak^ about two orders of mag- 
nitude slower than the main peak, has been hard to in- 
terpret by experimental means. Its strong dependence on 
the level of hydration, however, suggests that it should be 
linked to the protein. Indeed, our calculations give clear 
evidence that slow relaxation is related to overdamped 
motions of highly solvated ionized residues. The two or- 
ders of magnitude ratio of the slow and main relaxation 
times is in qualitative agreement with the dielectric mea- 
surements. Finally, the peak disappears when hydration 
of ionized residues is removed, which is analogous to dry- 
ing the sample in experiment. The assignment of the slow 
peak should also emphasize the involvement of water in 
the relaxation process. Since water is fast and follows 



FIG. 8. Superimposed structures of the bacterial enhancer- 
binding protein NtrC in dephosphorylated (light blue,PDB 
entry 1DC7) and phosphorylated (dark blue, PDB entry 
1DC8) states.—"^ The NMR structure was determined^ with 
carbamoylphospate binding to Asp54. The displacement of 
Glul24 in response to a probe charge at Asp54 is shown in 
Fig.H 

adiabatically the protein motions, ionized residues move 
by dragging hydration waters with them. 



C. Allosteric response 

The calculations of the allosteric response to ion bind- 
ing have been done for a single-domain signaling protein 
NtrC. The structures of this protein have been resolved^ 
both in unphosphorylated (denoted as NtrC) and in 
phosphorylated (denoted as P-NtrC) states. The latter 
state is short-lived, and it was maintained in solution at 
a large excess of the phosphordonor carbamoylphosphate 
carrying the charge of go 2 = —2. The addition of this 
charge to Asp54 active site (Fig. [8]) creates a Coulomb 
force acting on the neighboring atomic charges, such that 
each residue j experiences the force — (7o2Eoj. The per- 
turbing force produced by ion binding is therefore fairly 
nonlocal, in contrast to a common assumption,— and the 
calculation of the response requires full account of this 
fact. 

The frequency-dependent displacement of residue i in 
Eqs. ((26)) and (jST]) sums up all Coulomb forces acting on 
residues j from the active site labeled as "0". The un- 
phosphorylated (NtrC) state of the protein is very mo- 
bile, with several loops continuously changing their con- 
formation on the fis to ms time-scale. These motions 
mostly disappear in a more compact phosphorylated 
state.i^ Not surprisingly, we have found that sDENM 
is rather unstable and only DENM calculations could be 
done on the NtrC state. Therefore, sDENM calculations 
were done only on the P-NtrC state. The results for 
(5ri24(w) in both states are shown in Fig. [SI As expected, 
the inclusion of solvation in sDENM greatly enhances 
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FIG. 9. Frequency-dependent allosteric displacement [Eq. 
(pTll ] for residue i = 124 (Glu) of NtrC in the unphospho- 
rylated (NtrC) and phosphorylated (P-NtrC) states. The re- 
sults of calculations within DENM and sDENM are compared 
as shown in the plot. 




residue No. 

FIG. 10. Displacements Ar^ between two equilibrium struc- 
tures (NtrC and P-NtrC) of the NtrC protein. Sri{0) shows 
the displacement of residue i in response to placing a unitary 
probe charge = 1 at the position of C° of the binding site 
(Asp54, Fig.©. 



the displacement magnitude. The frequency dependence 
is also noteworthy. It implies the existence of /zs motions 
of the protein responding to the charge perturbation.— If 
the frequency of binding/unbinding events exceeds this 
frequency, the protein does not have the ability to re- 
spond to the perturbation and this subset of motions 
dynamically freezes. As a result, the displacement di- 
minishes. 

Figure [TO] emphasizes a strong effect of solvation on 
ion binding allostery presented in Fig. [9] by showing zero- 
frequency displacements of all residues in the NtrC pro- 
tein. The calculations have been done in DENM for NtrC 
and P-NtrC and in sDENM for P-NtrC. The results are 
compared with Ar^ displacements of C" between the two 
structures (Fig. [8]). It is clear that only by including sol- 
vation within sDENM does the displacements of residues 
in a given state reach the magnitudes comparable with 
the overall displacement amplitudes Ar;. It appears that, 
while solvation docs not strongly affect the global en- 
ergetics of ion binding (Fig. [3|), it critically affects the 
allosteric amplification of ion binding through conforma- 
tional transitions of individual residues. 



IV. SUMMARY 

Folded proteins have to maintain structural stability. 
At the same time many functions of enzymes and mo- 
tor proteins involve large-scale domain movements in re- 
sponse to binding and release of ligands. This require- 
ment makes one suggest that some stability needs to be 
sacrificed to allow amplification of a small perturbation 
into a large response. The question is what are the struc- 
tural motifs that allow amplification without compromis- 
ing the global stability. A related issue is the length of 
correlations, since long-ranged correlations are required 
for allosteric action at the distance. 

The first obvious target to address the problem is elas- 
ticity. The protein is densely packed and any force per- 



turbation propagates through its body as in a glass mate- 
rial. The elastic deformation spreads, however, through 
the elastic body and does not accommodate for a di- 
rected, specific action. Elasticity can capture motions of 
relatively rigid domains linked by flexible hinges )^ but to 
a lesser extent the allostery of monomeric, single-domain 
systems. 

One can alternatively turn attention to hydration 
water— as a possible medium for transferring signals. 
Water can store significant energy in the form of dipolar 
polarization and large entropy in its hydrogen-bond net- 
work, but it is also a highly non-specific medium. The 
required specificity might therefore reside at the protein- 
water interface combining large energies stored in hydra- 
tion with regulation achieved through identities of the 
surface residues. 

The surface charges, and to some extent dipoles, carry 
large solvation free energies and are strongly correlated in 
their motions, with water-mediated correlations decaying 
as r~^. Because of non-locality of correlations, an ensem- 
ble of ionized surface residues forms a strongly correlated 
net (or, in a sense, merged multiple pathways^) envelop- 
ing the entire protein. Given that multiple binding sites 
are typically involved in protein function, allostery might 
be designed not by building a specific site and attach- 
ing strings ("communication pathways"^^) to it, but by 
puling on the net wherever one finds a "knot". Some 
knots might be more important than the others from the 
perspective of biological function. The non-locality of 
this net excludes the possibility of well-defined pathways, 
they must be achieved by more specific interactions in- 
volving either no-polar residues^ or chains of hydrogen 
bonds 1^ 

A network of ionized, hydrated surface residues is a 
general property of all hydrated proteins, affecting a 
number of observable properties. The formalism of sol- 
vated dissipative electro-elastic network captures this re- 
ality and projects it on a number of susceptibilities de- 
scribing the response to a particular type of external per- 
turbation of a given experiment. A number of observ- 
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ables, such as rmsd's, dielectric response, electrostatic 
susceptibilities, and allosteric response are affected by 
solvation of the surface residues. The general outcome is 
that elastic motions of the residues become significantly 
heterogeneous, with softening achieved at the sites carry- 
ing charges. The interfacial heterogeneity does not dra- 
matically affect the global thermodynamics of the pro- 
tein or the thermodynamics of ion binding, but is criti- 
cal for local responses to external perturbations. While 
global motions of the protein altering its shape occur on 
the nanosecond time-scale, ^s motion o'^^'^^ are assigned 
to portions of the protein with highly hydrated ionized 
residues. 
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Appendix A: Solvation Integrals 

The solvation integrals in Eq. (|12|) are given by the 
following one-dimensional /c-spacc integrals involving the 
longitudinal structure factor S^{k) of the homogeneous 
solvent 

a{s,r) = Sr.o + - / Ji(sfc)2jo(rfc)(5^(fc) - l)dk, 

Jo 

b{s,r) = e{r - 2s) 

+ -(-) / Mskfj2{rk){S'^{k)-l)dk. 

(Al) 

Here, jn{x) is a spherical Bessel function, s denotes the 
effective radius of a residue, and r is the distance be- 
tween the centers of two beads in the network. Further, 
r = corresponds to one bead and that configuration is 
represented by the Kronecker delta Sr.o, which is equal 
to unity when r = 0. Since dipolar structure factors sat- 
isfy the asymptote S^{k) ^ 1 at fc — )• oo, this limit is 
separated from the numerical integral and is given by the 
first summand in each equation. 

The dipolar structure factors of polar liquids have been 
intensively studied in the past4^ Analytical models from 
liquid-state theories also exist4^ Several studies reported 
structure factors of force-field water modelsj ^^'^^ This 
function is accesible only from simulations since there 
is no known experimental technique giving access to it. 

For the purpose of estimating the integrals in Eq. 
(|A1[) we have taken the longitudinal structure factor 
of TIP3P water calculated from MD simulations;^ A 
simple parameterization of this function is available^^i^ 
based on the solution of mean-spherical closure for dipo- 
lar hard spheres4^ The results of this integration are 



shown in Fig. [TTJ As is seen the first integral involv- 
ing joirk) quickly goes to zero when reaching the dis- 
tance r/s ~ 2. Since r is either zero, for one-bead solva- 
tion, or greater than 2s, for different beads, only r = 
needs to be considered for this function. We therefore 
put a(s,r) = A{s)S^{0), where A{s = 4.4A) = 3.54 is 
numerically calculated. The value s = 4.4 A is the sum 
of the average radius of 3 A assigned to a residue and 1.4 
A for the radius of waters. 

The situation is just the opposite for the second inte- 
gral. It is zero at r = and reaches the value 5'''^(0) — 1 
at r/s = 2. This latter result implies that the continuum 
limit approximation S^{k) = S^{0) applies in this case. 
We therefore put 6(s,r) = S^{0)9{r ~ 2s). 

The overall result of these calculations, incorporat- 
ing solvent dipolar correlations through the longitudinal 
structure factor, is quite clear. The solvation energy, at 
r = 0, is renormalized by the factor A from the dielectric 
continuum limit S^{k) = S^{0). This renormalization 
effectively reduces the cavity radius for dipolar solvation 
from the distance of the closest approach of water to the 
residue s to s/A^^^. This is consistent with the com- 
mon observation that the effective cavity radius should 
fall between s and the van der Waals radius of the so- 
lute s — (Ts/S {<Js is the water diameter). On the other 
hand, the dipolar water-mediated coupling between dis- 
tant residues is well described by the continuum limit of 
the solvent dipolar response, and that fact is reflected in 
the constancy of 6(s, r). 
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